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ABSTRACT 

We study magnetic field steepening due to ambipolar diffusion (Brandenburg & Zweibel 1994) in 
protoplanetary disk environments and draw the following conclusions. Current sheets are generated in 
magnetically active regions of the disk where the ionization fraction is high enough for the magnetoro- 
tational instability to operate. In late stages of solar nebula evolution, the surface density is expected 
to have lowered and dust grains to have gravitationally settled to the midplane. If the local dust-to-gas 
mass ratio near the midplane is increased above cosmic abundances by factors > 10'^, current sheets reach 
high enough temperatures to melt millimeter-sized dust grains, and hence may provide the mechanism to 
form meteoritic chondrules. In addition, these current sheets possibly explain the near-infrared excesses 
observed in spectral energy distributions (SEDs) of young stellar objects. Direct imaging of protoplane- 
tary disks via a nulling interferometer or, in the future, a multi-band, adaptive optics coronagraph can 
test this hypothesis. 

Subject headings: instabilities — meteors, meteoroids — MHD — planetary systems: protoplanetary 
disks — solar system: formation — turbulence 



1. INTRODUCTION 

Chondrules are millimeter-sized silicate spheres with ig- 
neous textures embedded in the most common type of me- 
teorites, comprising 50-80 % of their mass. Their forma- 
tion is a long-standing mystery in solar system formation. 
Radioisotope datings place them among the first solids to 
have condensed out of the slowly cooling solar nebula. As 
such, they represent the remains of early materials that 
later formed the terrestrial planets. Their nearly spherical 
shapes and glassy compositions indicate that they cooled 
within a few hours after being molten, but not in the few 
seconds that objects that size would cool in free space. The 
associated heating event, whatever the origin, must have 
been transient and localized as evidenced by the inferred 
rapid heating and cooling rates, and the observed retention 
of volatiles such as sulfur (which vaporizes by prolonged 
heating above 650 K). To explain the high abundance of 
chondrules, it must have been a widespread phenomenon. 
Hence, encrypted in chondrules are important clues to 
physical conditions that prevailed in the early stages of 
solar system formation. Some even argue that taking the 
form of millimeter-sized molten spheres is a necessary in- 
termediate step for initially micron-sized dust grains to 
grow to kilometer-sized objects (Hewins 1997). 

A successful theory of chondrule formation should at 
least meet the following constraints: (1) the extremely 
short heating timescale of roughly seconds to minutes 
(Connolly & Love 1998), (2) the coohng rate of 10^-10^ 
K hr^^ (Radomsky & Hewins 1990; Lofgren & Lanier 
1990), and (3) the narrow range of sizes (0.1-1 mm), in- 
ferred from laboratory analyses of textures and composi- 
tions as well as from attempts to reproduce chondrules. 
In addition, many chondrules show evidences of multiple 



heating events, so the heating event must have occurred 
more than once. Since the initial discovery 200 years ago 
(Howard 1802), their formation process has received con- 
siderable attention in both meteoritics and solar system 
astronomy, yet it remains elusive (Wood 1996). Proposed 
mechanisms in the past (see the summary by Boss 1996, for 
the pros and cons) include nebula shocks, lightning, and 
an X-wind. None of these mechanisms have really been 
widely accepted (Hewins 1996), although recent work on 
shock heating by Desch & Connolly (2002) has gained wide 
attention. 

Infrared observations and subsequent theoretical mod- 
eling of T Tauri stars revealed the existence of circumstel- 
lar (protoplanetary) disks (Mendoza 1968; Shu, Adams, 
& Lizano 1987; Adams, Lada, & Shu 1987; Beckwith et 
al. 1990), many of which are actively accreting (see, e.g. 
the review by Calvet, Hartmann, & Strom 2000). These 
systems are generally thought to be analogs of the early 
solar nebula. On the theoretical side, magnetorotational 
instability (MRI) provided a solution to the long-standing 
problem of identifying the cause of anomalous viscosity, 
which allows mass accretion of the disk (Balbus & Hawley 
1991; Hawley & Balbus 1991, collectively referred to as 
BH91 hereafter). When applied to partially ionized disks, 
it remains an efficient mechanism of generating turbulence 
as long as ions are sufficiently well-coupled to neutrals, i.e., 
if ion-neutral collision rate far exceeds the orbital rotation 
rate: t'in ^ fl (Blacs & Balbus 1994; Mac Low et al. 
1995; Hawley & Stone 1998). In the case of protoplane- 
tary disks, however, a large portion, especially near the 
disk midplane, is associated with relatively low ionization 
fraction (~ 10~^^), and hence thought to be magnetically 
dead (Gammie 1996; Glassgold, Najita, & Igea 1997). 



In partially ionized disks, ambipolar diffusion arises be- 

^ Department of Astronomy, Columbia University, 550 West 120th Street, New York, NY 10027; moo@astro.coIumbia.edu 

^ Department of Astrophysics, American Museum of Natural History, 79th Street at Central Park West, New York, NY 10024; morde- 
cai@amnh.org 

^ Department of Earth and Planetary Sciences, American Museum of Natural History, 79th Street at Central Park West, New York, NY 10024; 
debel@amnh.org 



2 



Joung, Mac Low, & Ebel 



cause magnetic forces act on the charged component, i.e. 
ions, only (Mestel & Spitzer 1956). As the ions drift 
through neutrals, the magnetic field lines drag the ions 
with them, steepening field gradients near magnetic nulls 
stretched by shear flows (Brandenburg & Zweibel 1994, 
hereafter BZ94). Consequently, sheets of high electric cur- 
rent form, as was confirmed by numerical simulations (Mac 
Low et al. 1995; Suzuki & Sakai 1996). BZ94 originally 
proposed this as a mechanism for rapid magnetic rccon- 
nection in the interstellar medium. Ambipolar diffusion is 
important in another context: as a heating source (Scalo 
1977; Goldsmith & Langer 1978; Zweibel & Josafatsson 
1983; Elmegreen 1985). Recently, using three-dimensional 
simulations of turbulent, magnetized molecular clouds, 
Padoan, Zweibel, & Nordlund (2000) pointed out that am- 
bipolar diffusion heating is a significant source of heating 
in molecular clouds and that it can be larger than the cos- 
mic ray heating. In this paper, we examine the possibility 
that similar physics operates in protoplanetary disks, such 
as in the solar nebula that eventually became our solar 
system. 

The purpose of this paper is twofold. First, we show that 
current sheets are generated in magnetically active regions 
within protoplanetary disks. We study the properties of 
these current sheets and determine when and where the 
associated heating becomes significant. Second, we apply 
this formalism to the following two unresolved issues: the 
heating mechanism that melted chondrules in the early so- 
lar nebula and the near-infrared excesses observed in spec- 
tral energy distributions (SEDs) of yoimg stellar objects. 

We later demonstrate that the heating of dust grains by 
current sheets is ineffective for a standard model of the 
early solar nebula in which the dust-to-gas ratio C — 0.01. 
However, at lower gas densities and much higher dust- 
to-gas ratios than those usually assumed in the standard 
model, these current sheets are found to reach high enough 
temperatures to melt millimeter-sized dust grains. We 
propose that these unusual conditions occur near the mid- 
plane in late stages of solar nebula evolution, when the 
gas surface density has decreased due to various mecha- 
nisms discussed below, and dust grains have gravitation- 
ally settled. The surfaces of this dust-enriched midplane 
layer are then sufficiently dense and ionized to host thin 
{I « 10^ km), hot (T > 1600 K) current sheets. Circula- 
tion of dust to the surface by turbulence within the dust 
layer is invoked to explain the observed high abundance of 
chondrules. 

The basic features of our current sheet model are set 
forth in § 2. In § 3.1 we use this model to investigate 
a scenario in which chondrules are heated inside current 
sheets, solving numerically the coupled equations of en- 
ergy transfer for gas and dust particles. The results are 
presented in § 3.2. We argue, in § 3.3, that in late stages 
of solar nebula evolution current sheets capable of melt- 
ing most of the millimeter-sized dust particles form in the 
dust-rich midplane. In § 3.4 we discuss possible observa- 
tional signatures of current sheets. Finally, in § 4 we list 
limitations of our model. 



2. CURRENT SHEET MODEL 

2.1. Basic Equations 



The magnetic field evolves according to the induction 
equation. 
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where v,, denotes the electron velocity and is the associ- 
ated resistivity given by 

_ (? _ 234 

47r£7e Xe 

(Krall & Trivelpicce 1973). Here, c is the speed of light, 
Xe = n,./nn is the electron fraction, and Tg is the ambi- 
ent gas temperature. For the last equality in Eqn. (2), 
we adopt the value of the electrical conductivity, a^, from 
Draine, Roberge, & Dalgarno (1983) and Blaes & Balbus 
(1994). Subscripts e, i, and n denote quantities pertain- 
ing to electron, ion, and neutral components, respectively. 
Other quantities of relevance in this paper are the mean 
mass per particle of m„ = 2.33 rup, where rrip is the proton 
mass, and the mean ion mass of rrii = 39 mp at densities 
where ionized potassium, K+, is generally the dominant 
ion species (Balbus & Tcrquem 2001). 

In the ideal MHD limit, the electron velocity is indis- 
tinguishable from the center-of-mass velocity of the fluid. 
This approximation breaks down in two ways as Xe de- 
creases. First, ions slip through neutrals, and for low 
enough Xei when ion pressure and inertia can be neglected, 
the strong coupling approximation can be used to com- 
pute the ambipolar drift velocity by balancing the Lorentz 
force J xB/c acting on the ions against the frictional drag 
{vi - Vn) arising from coUisions with neutrals. Sec- 
ond, when ions and electrons move along a direction per- 
pendicular to a magnetic field, they experience the Lorentz 
force acting normal to both directions but they drift away 
in two opposite directions due to the different signs of their 
electric charges. This causes the Hall effect. The electron 
velocity in a weakly ionized medium is then 

Ve = V + {Vi - V) + {Ve - Vi) = V + ^ ^ ^ ^ ^ (3) 

where pi = mini and Vin = "fPn are the ion mass density 
and ion-neutral collision rate. The value of the drag coef- 
ficient, 7 = 2.75 X 10^^ cm^ s~^ g~^ (Balbus & Terquem 
2001). Since Xe ^ 1 in our case (see § 2.2), the neutral 
velocity is the same as the center-of-mass velocity, v. The 
cmrent density has the usual form J = (c/47t)V x B. We 
plug Eqn. (3) into Eqn. (1) and consider the case in which 
the currents and fields are perpendicular (J • B = 0), so 
the induction equation becomes 
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The Hall term does not affect the time evolution of B, 
since Vx(JxB) a — Vx VB^ vanishes. Thus, ambipolar 
diffusion acts as an additional field- dependent resistivity. 
We define the ambipolar diffusion coefficient as 

>^AD = (5) 
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where Bq means B in gauss. 

For our model of the protoplanetary disk, we take 
the minimum- mass solar nebula of Hayashi (1981). It 
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is cliaractcrizccl by the siirfacc mass density S = 
Eo (i?/3 AU)-3/2 with Sq = 3.3 x 10^ g cm-^, where R 
is the radial distance from the central star, and the tem- 
perature profile given by Tg = 1.6 x 10^ (i?/3 AU)"^/^. We 
expect So to decrease in time as the nebular gas clears out 
in the late stage of protoplanetary disk evolution. Hence 
we use the form Sq = 3.3 x lO^/s g cm"^, where /e is 
an arbitrary parameter less than or equal to unity. If we 
take the gas to be isothermal, the vertical density profile 
in hydrostatic equilibrium has a Gaussian form (Kenyon 
& Hartmann 1987; Chiang & Goldreich 1997, hereafter 
CG97) 



no exp 



(6) 



where no = E/(\^m„/i) = 1.4 x 10^^ h (ii/3AU)-iV4 

cm~^ is the number density of gas in the disk midplane, z is 
the vertical distance measured from the midplane, and h = 
Cs/fl is the scale height. Here, as usual, Cs = {^Tg/nin)'^^'^ 
is the sound speed and f2 is the Keplerian rotation rate. 
With these values, we can roughly estimate the impor- 
tance of ambipolar diffusion by taking the ratio of the two 
coefficients in Eqn. (4): 
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We note that (1) with a magnetic field sufficiently ampli- 
fied (up to -^1 G), ambipolar diffusion is at least as im- 
portant as Ohmic dissipation and (2) the ratio has strong 
spatial dependence within the disk. Remnant magnetism 
in meteorites siiggcsts that _B w 1 G was typical in the 
solar nebula. If an efficient dynamo is in action, the am- 
bipolar diffusion term would generally be important in the 
outer (i? ^ 1 AU) region of the disk. In the inner region, 
for a similar value of B and /s = 1, it is nonnegligible 
only a few scale heights above and below the midplane, 
where gas density is much lower. It generally becomes 
increasingly important as /e decreases. 

2.2. Magnetorotational Instability in Protoplanetary 

Disks 

Compared to other astrophysical systems, protoplane- 
tary disks are unusually cold and dense and are there- 
fore associated with very low ionization fractions. Ioniza- 
tion fraction is a key parameter in determining whether 
a disk will be turbulent or not, since the MRI is confined 
to regions where the ionization fraction is high enough 
to ensure good coupling between ions and neutrals. This 
occurs when the ion- neutral collision rate (vm) is signifi- 
cantly greater (by 3 100 at least) than the rotation rate 
(n) (Blaes & BallDus 1994; Mac Low et al. 1995; Hawley & 
Stone 1998). The critical ionization fraction is calculated 
to be 

(e.g., Balbus & Hawley 1998). Interestingly, this value is 
only slightly higher than the values of 10~^2-10~^^ esti- 
mated in the bulk of protoplanetary disks (Gammie 1996; 
Glassgold, Najita, & Igea 1997; Sano et al. 2000). 



Many numerical simulations have demonstrated that 
MRI-induced velocity shears amplify the magnetic field 
(see, e.g., Hawley & Balbus 1991; Brandenburg et al. 1995; 
Hawley, Gammie, & Balbus 1995). In our model, per- 
turbations in the magnetic field are also assumed to be 
generated by the MRI. Considering the low values of Xi 
estimated in protoplanetary disks, it is initially unclear 
whether this assumption is valid or not. Therefore we ask: 
which region of the solar nebula, if any, is prone to the 
MRI? (Previous numerical simulations assumed ionization 
fraction Xi ^ Xcrit-) This question has been investigated 
by several authors (Umebayashi 1983; Sano et al. 2000; 
Fromang, Terquem, & Balbus 2002). Their results are 
as follows. The ionization fraction is sufficiently high in 
the innermost region of the disk up to a few tenths of an 
AU, where the ambient gas temperature is high enough 
{Tg > 10^ K) for coUisional ionization to be effective 
(where it dominates over cosmic ray and radioactive el- 
ement ionizations; Umebayashi 1983). Also, cosmic rays 
and stellar X-ray radiation ionize the surface layers (with 
thickness given by the stopping depth of Sa w 10^ g/cm^) 
at R<10 AU and the entire disk at i? > 10 AU, where the 
surface density of the disk falls below S^. These, then, are 
the regions where the disk is turbulent, hence effectively 
transferring angular momentum and also actively accret- 
ing. On the other hand, the region near the disk midplane, 
where terrestrial planets and the asteroid belt now lie, is 
generally thought to contain no significant heating source 
and be quiescent, i.e. magnetically dead (radioactive heat- 
ing is too inefficient to maintain Xi > Xcrit] see Gammie 

1996) . These considerations led to the layered accretion 
model of Gammie (1996; see also Glassgold, Najita, & Igea 

1997) . 

There are two length scales, other than the scale height 

h, of relevance to this discussion. The wavelength of per- 
turbation most unstable to the MRI is only slightly larger 
than the critical wavelength at which the instability sets 
in (BH91) and is given by 



= -^«0.23/,-.. 



3Auj 
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AU 
(9) 

(BH91). Other processes affect the growth of perturba- 
tions: in dense parts of the disk, collisions with neutrals 
hinder the charged components from drifting with the 
magnetic field. In this case, Ohmic dissipation acts to 
stabilize wavelengths shorter than the resistive diffusion 
scale An = 2-Krj/vA, where va = B/{A'KpnYI'^ is the Alfven 
speed. In other words, the critical wavelength at which 
MRI sets in is given by \crit = fnax {Xbh , ^n) (Blaes & 
Balbus 1994; Sano et al. 2000). When Xn > Xbh, the 
perturbation growth becomes much slower. An additional 
constraint enters when Xcrit becomes comparable to the 
scale height h; at this point, satmation of the MRI is 
thought to occur (Balbus & Hawley 1998). 

Sano et al. (2000) investigated where the MRI operates 
in protoplanetary disks by including a detailed chemical 
reaction scheme and the effect of the recombination of 
ions and electrons on grain surfaces. They showed that 
for the minimum-mass solar nebula (/e = 1) the 2-3 AU 
region of the disk would be stabilized by Ohmic dissipa- 
tion. However, as the surface density parameter, /e, is 
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lowered, the unstable region expands (see their Figure 8). 
Also, as the disk evolves, dust grains settle to a thin layer 
in the midplane, and so the grain abundance in the bulk 
of the nebula gradually decreases. For this reason, Sano 
et al. (2000) introduced a dust depletion factor, fd, which 
measures the abundance of dust grains relative to that 
in molecular clouds. Thus a low value of fa may corre- 
spond to an old, evolved disk. They demonstrated that 
as fd decreases for a given surface density, the stable re- 
gion (the "dead zone") shrank in size, and for fd < 10~^, 
the entire disk becomes unstable for i? > 3 AU (see their 
Figure 11). This happens because the recombination rate 
on grain surfaces decreases. They found a similar trend as 
grain growth occurs. From these results, we conclude that 
the region at i? « 3 AU will become magnetically active 
at a sufficiently late stage of the evolution. We use and 
discuss this result further in § 3. Figure 1 compares the 
three relevant length scales at an early (a) and a late stage 
(b) of solar nebula evolution. The region of interest where 
magnetic fields amplify due to the MRI extends from 2 
AU to 4 AU in (a), while with a reduced gas density, the 
active region moves inward to 1-2 AU (b). 

2.3. Formation of Current Sheets 

In this and the following sections, we consider a simple, 
one-dimensional geometry illustrated in Figure 2. Also 
shown are the directions of various terms in Eqns. (3) 
and (4). Consider first a case in which a magnetic null 
is present. Near the null, a magnetic pressure gradient, 
— V(i3^/87r), acts on ions and pushes the field lines into 
the null point in both directions. The ion-neutral drift is 
then directed toward the null. Along a gradient in the 
magnetic field, this nonlinear term produces faster field 
diffusion in regions of stronger field, driving field into re- 
gions of weaker field and steepening the gradient, rather 
than spreading it out as the linear Ohmic term docs. 
Ultimately, the gradient sharpens towards a singularity, 
marked by a sheet of high electric current. This steepening 
of magnetic field profile (BZ94) occurs within a timescale 
of TAD = -^^B///(2 Aad) = (ti"^/ 6) (j/i„/0^), calculated to 
be 7r/12 of the orbital period when Xi = Xcrit- The factor 
of 1/2 in TAD takes into account the fact that the charac- 
teristic length decreases as steepening proceeds. 

Remarkably, this sharpening has been shown to occur 
even in the absence of nulls due to a variety of processes 
including a rotating eddy (BZ94; Zweibel &: Brandenburg 
1997), accretion disks sustaining the MRI (Mac Low et 
al. 1995), large amplitude Alfven waves (Suzuki & Sakai 
1996), and the nonlinear stage of Wardlc C-shocks (Mac 
Low & Smith 1997). As previously mentioned, we take 
the example of magnetohydrodynamic turbulence driven 
by the MRI. In this case, velocity shears produce magnetic 
nulls and current singularities . Stretched by shear flows, 
magnetic nulls develop where the magnetic field changes 
directions. We therefore crudely estimate that the separa- 
tion between two adjacent magnetic nulls Lbh = ^bh/^- 

We note that ambipolar diffusion allows the field to slip 
through the neutrals and tends to weaken the MRI. For 
strong enough diffusion, i.e. when the ion-neutral colli- 
sion rate {vm) is much higher than the growth rate of 
the MRI (~ f2), the instability is entirely prevented (Mac 
Low et al. 1995; Hawley & Stone 1998). On the other 



hand, ambipolar diffusion is required to generate current 
sheets. For our model, we use parameters on the inter- 
face between the two regimes. This choice is justified by 

two-dimensional simulations performed by Mac Low et al. 
(1995), which showed that for drag coefficients 7 a few 
times the minimum value that allows the instability to de- 
velop, jmin = ^ / {mi Un Xcrit), both the instability and 
ambipolar diffusion sharpening occur. In other words, we 
require 1 < "//"/min < 100. Alternatively, magnetic field 
amplification could have preceeded the epoch of current 
sheet formation. 

We now make quantitative estimates of the properties 
of current sheets. BZ94 showed that, in the absence of 
resistivity, the magnetic field profile relaxes to i? oc x^^^, 
and the current density J oc dB/dx oc x~'^/^ becomes 
singular at the origin. In reality, both the finite resis- 
tivity and the ion pressure act to remove this singularity 
(Brandenburg & Zweibel 1995). Under the physical con- 
ditions of protoplanetary disks, ion pressure can be ne- 
glected (Heitsch & Zweibel 2003). The induction term in 
Eqn. (4), V x {v x B), becomes progressively less im- 
portant as the field steepens; I/O ~ vL/r] where / and 
O denote the magnitudes of the induction and Ohmic dis- 
sipation term respectively, and is therefore ignored. The 
quasi-steady-state magnetic field profile ( "Quasi" since the 
magnetic field profile lasts only as long as the velocity 
shear is present, for about a dynamical time) is calculated 
by balancing the Ohmic dissipation term against the am- 
bipolar diffusion term in Eqn. (4): 

(S» + «J)fg. + 2B.(^y = 0. (10) 

Here, for simplicity, we use the definition 

Bo = (47r77p,i.„)i/2 « 0.26 ngi3 t'J^ G , (11) 
which measures the relative importance of Ohmic dissipa- 
tion to ambipolar diffusion. We used n^is and Tgs to de- 
note ng/(10^^ cm~^) and Tg/(10^ K), respectively. Ohmic 
dissipation dominates where B < Bq] ambipolar diffusion 
wins where B > Bq (sec Eqn. 7). Remarkably, Bq is inde- 
pendent of the ionization fraction; so is the magnetic field 
profile. Solving Eqn. (10), we obtain 

^B^y+B'oBy = rx, (12) 

assuming that the field vanishes at the origin. An in- 
tegration constant, F, is determined by the boundary 

condition By(x = Lbh) = Bmax, which leads to F = 
BmaxiBo + B'^^^/3)/Lbh- The ampfitude of the mag- 
netic field, Bmax, is determined from the remnant mag- 
netic field strengths recorded in the matrix of chondrites 
(0.1-1 G) and in chondrules themselves (1-10 G) (Desch & 
Connolly 2002, and references therein). From Eqn. (12), 
we see that Bq ^ near the origin, and the effect of 

resistivity yields By « Fx/Bq. Note the current density 
singularity is now removed; the maximum current density 
near the origin is given by J{x = 0) = cT /AttB^, which is 
finite. Sufficiently far from the origin where Bq <C B^/3 is 
satisfied. By « (3ra;)^/^ holds, and we retain the result in 
the non-resistive limit of BZ94. There, the current density 
declines steeply (J oc x~^/^). We therefore set B^ = By/3 
and compute the width of the current sheet to be 

U = ^ (13) 
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This expression is valid for Bmax > V^Bq, which is satis- 
fied under the conditions of interest. 

2.4. Heating Mechanisms 

Inside current sheets, two heating mechanisms act. 
First, electrons and ions arc accelerated by the electric 
field present, and then collide with neutrals and heat up 
the entire gas, leading to familiar Ohmic dissipation. The 
corresponding heating rate per unit volume essentially fol- 
lows from the magnetic field profile: 



j^^±(dBy 
(Te 47r \ dx 



(14) 



There is additional heating due to friction when ions, accel- 
erated by the Lorentz force, diffuse through neutrals (see, 
e.g., § 2.1 of this paper; Padoan, Zwcibel, & Nordlund 
2000). In this case of ambipolar diffusion heating, 



Sad = Pil^in\Vi - Vn\ = 



nBl 



rjBl [dB. 



dx 



(15) 



yields the volumetric heating rate. Note that their ratio 
is given simply by EadI^q. = By/4,BQ. Within the range 
of parameters we are interested in, the ambipolar diffu- 
sion heating integrated through the current sheet along 
X exceeds (by a factor of ~10 in our fiducial model) the 
contribution from Ohmic dissipation. 

The resulting heating profiles (log e vs. log x) as a func- 
tion of Bf), defined in Eqn. (11), arc displayed in Figure 
3. They are naturally symmetric around x = Q. In terms 
of shape, these double peaks are similar to what has been 
observed in numerical simulations of Padoan, Zweibel, & 
Nordlund (2000), although their sheets are produced by 
supersonic shocks and density compressions. Ambipolar 
diffusion sharpening becomes more efficient at lower gas 
densities (sec Eqn. 13), and the current densities become 
close to singular. This implies stronger and more concen- 
trated heating far above and below the disk midplanc. Al- 
ternatively, if the disk becomes thinner with time (Ruden 
1999), the heating rate also increases. For example, soon 
after thermonuclear reactions began in the Sun, a power- 
ful solar wind is thought to have blown away the gas ~0.1 
to a few million years after the solar nebula formed. Ac- 
cretion onto the central star or, in local regions, planetary 
accretion as gap formation occurs would have the same 
effect. A density drop caused by any of these mechanisms 
introduces a regime where chondrule formation is favored. 
As the overall gas density Ug drops in time, Bq oc Ug, so 
the peak heating rate rises. 

3. APPLICATIONS 
3.1. Energy Equations and Chondrule Formation 

We investigate chondrule heating in current sheets by 
solving numerically the coupled equations of energy trans- 
fer for gas and dust particles, and compare the results with 
the experimental constraints on chondrule formation. 

Dust particles are usually assumed to possess a power- 
law distribution of sizes with upper and lower size cutoffs 
(e.g. Chiang et al. 2001). For simplicity, we use just two 
types of dust particles: small, micron-sized dust grains 



with radius Od — 1 /^m and bigger, millimeter-sized chon- 
drules with Oc = 1mm. Also, it is assumed that the 
same amount of mass is contained in each type. Near the 
dust melting temperature, small grains radiate in the near- 
infrared most efficiently. Hence their thermal evolution is 
of considerable importance to the overall development of 
current sheets. 

To study how these current sheets affect the tempera- 
tures of dust particles that emit the bulk of the infrared 
radiation, we use the one-dimensional slab geometry con- 
sidered in Hood & Horanyi (1993). Dust and gas travel 
through a current sheet with relative velocity of Vr, where 
Vr is a free parameter of the model but assumed to be 
comparable to the Alfven speed, va = i3/(47rp„)^/^, and 
smaller than 17 km s~^, the Keplerian orbital speed at 
R = 3 AU. For our fiducial model, Vr = 1 km s~^. We 
discrctizc the region for computation so that each zone has 



lateral extent of our computational domain exceeds that 
of our fiducial current sheet {Ics ~ 10^ km). The gas and 
dust temperatures are uniform in each zone but in general 
differ from each other. The temperatures evolve in time 
and are assumed to vary only in one spatial direction (x) . 

Entering a current sheet, gas molecules gain energy 
from the two heating mechanisms mentioned in § 2.4 
and lose it through collisions with dust grains and chon- 
drulcs. The coUisional heat transfer rate per unit surface 
area between gas and dust, denoted by q^, is defined as 
qd = (3/4-v/7r) {2k/mnf^'^ Pg{Tg-Td) and qc is defined sim- 
ilarly for chondrules (Hood & Horanyi 1993). Then the gas 
temperature evolves in time according to 

dT„ 



dt 



■ ^~ m '^^"■i 1i+P9^iR i-^r - 4:CrT^) , (16) 



where k is the Boltzmann constant, rij is the number den- 
sity of either dust grains (i = d) or chondrules (i = c), and 
e = en + ead- The opacity of the gas, k/_r, is mainly due 
to dust associated with it, and is set to 1.0 cm^ g""'^, appro- 
priate for infrared wavelengths (Chiang et al. 2001). The 
intensity of radiation, Jr (integrated over all wavelengths) 
at x = X{) from plane-parallel, temperature-stratified at- 
mospheres has the following form (Mihalas 1978; Desch & 
Connolly 2002): 

r+oc 

Jr{xo) =2a {pgKiRT^ 



Hi naf €iT*) Ei{P\x—xo\) dx 
(17) 



with /? = PgKiR Tg + J2i=d c '^i ""'^i ^i- Here, Ci is the emis- 
sivity of a grain of type i averaged over the Planck function 
at temperature Ti{x), which wc simply set at 1.0 (for com- 
parison, Desch & Connolly (2002) uses values of 0.8 and 
0.4 for 1 mm- and 1 /im-sized grains, respectively) , and Ei 
is the exponential integral of the first type. Dust grains 
and chondrules, on the other hand, arc heated both by col- 
liding with hotter gas molecules and by absorbing radia- 
tion from nearby smaller, micron-sized, grains which reach 
high temperatures quickly. They lose energy by radiating 
an approximate blackbody spectrum at Td and Tc, with 
emission efficiency of Cd and Cc for dust and chondrules, 
respectively. Therefore, 

dT- 

miCi-^= Tra^ [Aqi + ei ( J. - AaT^)] • (18) 
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As before, the index i can be either d (for dust) or c (for 
chondrule). The last term in Eqn. (18) represents the in- 
cident flux on a given particle radiated by dust grains in 
the other zones minus its own emission. 

The experimental constraints on chondrule formation 
are (1) the grains are completely melted for a few minutes; 
(2) the cooling timescale of gas and dust is approximately 
an hour, consistent with the oft-quoted value of 10^-10^ K 
hr^^; and (3) small grains (< 0.1 mm) are absent, possibly 
evaporated or destroyed completely after passing through 
a current sheet, leading to the observed narrow size distri- 
bution centered at ~ I mm. 

From the properties of current sheets described in § 2.3- 
2.4, we infer that the above constraints on chondrule for- 
mation will be satisfied only if the following conditions are 
met: 

(1) In order to obtain high enough heating rate, the gas 
density has to be low {ug < 10^^ cm^''), because the heat- 
ing rate increases steeply as Ug decreases. 

(2) To have "localized" heating events, spatially averaged 
chondrule/ dust densities should be high (n^ > 1 cm~^). 
Otherwise, the mean free path for IR photons will exceed 
a few hundred kilometers and the heating timescale will 
be much more than "a few minutes." 

(3) To ensure good thermal coupling between gas and dust, 
either the gas density has to be high {ug > 10^^ cm~^) be- 
cause Qd, Qc oc Hg which cannot be satisfied because of 
the first condition, or dust grains need to be very small 
(ttd < 1 /um) so they approach Tg quickly. 

In order to satisfy conditions (l)-(3) simultaneously, the 
dust-to-gas mass ratio, C, should be very large, about 50. 
Usually C « 0.01 is assumed for the minimum mass solar 
nebula, but after a few million years the gas density is ex- 
pected to have decreased due to viscous diffusion (Ruden & 
Lin 1986) as well as the various mechanisms mentioned in 
the last paragraph of § 2.4. Dust settles into a thin, dense 
layer on a timescale sa (a/1 ij,m)~^f^ Myr, independent 
of R (Ruden 1999, CG97), comparable to the lifetimes of 
protoplanetary disks. There will still be some dust grains, 
particularly submicron-sized ones, that remain well above 
the midplane mixed with the gas. But at late times the 
dust density outside the thin, dense dust layer is likely to 
be quite small, decreasing the dust-to-gas ratio there. In- 
side the dust layer, however, solids may dominate the local 
density (Cuzzi, Dobrovolskis, & Champney 1993). Hence, 

« 50 may be reached in the midplane, where most solids 
reside. 

We solve Eqns. (16)-(18) simultaneously for Tg{x,t), 
Td{x,t), and Tc{x,t) using a fifth-order Runge-Kutta 
scheme appropriate for a stiff set of equations with adap- 
tive stepsize control (Press et al. 1992). Both gas and dust 
particles begin at the same temperature, Tj„jt = 500 K. 
The center of a current sheet is initially placed at a; = 1000 
km and moves to the left at Vr- To simplify the problem, 
we ignore the time evolution of e and Xi (see § 4). For our 
fiducial model, rig = 10^^ cm""^ and both and ric are 
set to yield desired values of the dust-to-gas ratio C- 

3.2. Results 

Figure 4 shows the time evolution of gas, dust, and chon- 
drule temperatures at a fixed position as a current sheet 
moves through, for our model runs with ( = (5, 50). We 



also ran the code for ( = 0.05 but do not plot the result 
because both T4 and Tc remain practically unchanged at 
Tinit for over 10^ s. We find that for ^ = 5 the gas temper- 
ature rises appreciably (above 1000 K) but, due to poor 
thermal coupling between gas and dust, neither nor Tg 
changes significantly, reaching only ~ 700 K after 10^ s. In 
the run with 1^ = 50, however, due to the increase in optical 
depth, dust grains follow the gas temperature closely and 
the radiation field also steadily increases in strength. The 
grains melt after about 570 seconds and evaporate after 
about 630 seconds, at which point we stop the simulation, 
as it does not take into account the evaporation of par- 
ticles. We changed the boundary condition between fiat 
at Tinit and a floating boundary condition but found the 
result remained qualitatively the same. In an improved 
model that accounts for evaporation and a range of dust 
grain sizes, wc expect the smallest grains to evaporate first, 
decreasing the optical depth, thus widening the heated re- 
gion and lowering the net volumetric heating rate. The 
larger grains will therefore probably not evaporate. This 
is interesting since it could naturally explain the peaking 
of the size distribution of chondrulcs around 1 mm. 

The fact that chondrules get melted only for high val- 
ues of C is consistent with the FeO record in chondrules, 
which seems to indicate oxygen fugacities (relative abun- 
dance of oxygen to hydrogen) at formation that are higher 
than solar values by at least two or three orders of magni- 
tude (Wood 1967). This is most often interpreted to mean 
that chondrules were heated in regions with extremely high 
dust-to-gas ratios. The dust is taken to be oxides of Si, 
Mg, Ca, Al, etc. in chondritic proportions which bring 
in a large amount of oxygen. For this reason also, it is 
advantageous to evaporate the smaller grains first because 
then oxides will be incorporated into the gas and the par- 
tial pressures of condensable elements will increase locally. 
This helps to make liquids stable at high temperature, and 
also suppresses the evaporation of silicate liquids which 
would have caused heavy isotope enrichment of resulting 
chondrules; such enrichment is not observed (Galy, Young, 
Ash, & O'Nions 2000). Ebel & Grossman (2000) gives a 
full treatment of what dust enrichments are necessary to 
make FeO-bearing silicates and also to make chondrule liq- 
uids stable at low pressure relative to solids, as well as a 
discussion of dust-enrichment scenarios with high oxygen 
fugacities. The idea that the midplane is a better place 
to consider chondrule formation is supported again in this 
context. 

The value of the gas density used in Figure 4 is = 10"'^^ 
cm^"^ corresponding to /s = 10^^. Since we expect the 
gas density in the disk to decrease to this level after a time 
of order a few million years, this argues for chondrule for- 
mation in a late, evolved disk, and is consistent with the 
age difference between CAIs (refractory Ca-Al-rich inclu- 
sions) and chondrules, estimated from the decay of ^^Al 
(with a half-life of 0.7 Myr) to be 2-3 Myr (Hutcheon et 
al. 1994; Russefi 1996; Swindle et al. 1996). 

3.3. Current Sheets in the Disk Midplane 

From the results in the previous section, we conclude 
that in a standard model of the early solar nebula in which 
dust and gas are well-mixed, the disk will have ( « 0.01 
everywhere and hence heating of dust grains by current 
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ing shows that the turbulence diffuses inward and per- 
sists even very near the midplane (Cuzzi, Dobrovolskis, 
& Champncy 1993, and references therein)^. We take a 
typical velocity of U ^ 10^ cm s^^ [sec Eqn. (65) of the 
previous reference] and the thickness of the dust layer, 
Ld PS 3 X 10'' km, and estimate the eddy turnover time 
to be Ld/U « 1 yr. Keeping in mind that the fraction 
of current sheets covering the dust surface is quite small 
{Ics/Lbh ~ 3 X 10^"''; sec the next subsection), we expect 
an average dust grain to be exposed to current sheet heat- 
ing about once every 3 x 10* yrs. After the brief heating, 
now spherical, mm-sized chondrulcs will settle down to 
the midplane, collecting rims of dust in the process. Such 
rims are observed (Metzler & Bischoff 1996, and references 
therein) . 

3.4. Near- Infrared Excesses 

The higher than predicted flux at 2-8 /xm observed in 
the spectra of some young protoplanetary disks has been 
known for a decade (Hartmann, Kenyon, & Calvet 1993; 
Natta et al. 1999; Chiang ct al. 2001). Using a disk model 
that included an accounting of grain size distributions and 
grain compositions, and numerical solution of the equa- 
tions of radiative and hydrostatic equilibrium under the 
two-layer approximation of CG97, Chiang et al. (2001) 
rediscovered this in 4 of their 5 Herbig Ae and T Tauri 
stars, and found a relative dearth of emission shortward 
of 10 /xm arising partly from low absorption efficiencies of 
silicate grains. It is noteworthy that these sources also 
showed evidence for dust settling; parameters appropriate 
for gas and dust that arc well-mixed in interstellar pro- 
portions did not fit the overall level of infrared excess at 
A < 100 nm for 4 of them. They interpreted this to mean 
that dust in disk surface layers has settled vertically to- 
wards the midplane. 

DuUemond, Dominik, & Natta (2001) proposed a disk 
model in which an inner rim located at the dust evapora- 
tion radius forms a vertical wall due to direct heating by 
the central star and. as a result, puffs up and intersects 
more starlight. Its reemitted flux produces a bump that 
peaks at 2-3 /zm in the SED. They note that this modifi- 
cation would be significant only for Herbig Ae stars. This 
mechanism works only if optical depth due to gas opacity 
in the inner hole is negligible, but in fact it is quite large, 
of order 10^, when a standard surface mass density profile 
S oc is assumed, although the surface density of 

the inner disk is not constrained by the data (Natta et al. 
2001). 

We suggest an alternative mechanism by noting that 
current sheets at the dust melting temperature would radi- 
ate at just the right wavelengths to fill in the deficient flux 
at 2 -8 yum, perhaps explaining the excess near-IR fluxes. 
A simple estimate can be made as follows. We take the 
mean separation between current sheets (produced by the 
Brandenburg- Zweibel mechanism) to be Lbh = ^bh/'2', 
roughly the distance between two adjacent magnetic nulls 
produced by the MRI. The total emission from dust grains 

* The dominant ion species in protoplanetary disks are alkali metals (Na+, K+). Grains are dominant charge carriers at Tg < 103 K, but at 
higher Tg desorption of charged ions from the surface of dust grains becomes effective (Umebayashi 1983). Thus, our neglect of the effect of 
charged dust grains is partially justified, at least inside current sheets. 

® But see Sekiya (1998) and Youdin & Shu (2002) for claims that turbulent mixing becomes ineffective if the dust-to-gas surface mass ratio is 
significantly enhanced. 



sheets will not be effective. This conclusion no longer holds 
in the late stage of the disk evolution. As argued by Sano 
et al. (2000), the bulk of the disk becomes unstable to 
the MRI, therefore susceptible to current sheet formation, 
for low enough values of the surface density factor /s or 
of the dust depletion factor because the recombination 
rate on grain surfaces decreases and thus Xi rises except in 
a thin layer around the midplane. Current sheets capable 
of melting mm-sized dust grains form only when rig < 10^^ 
cm~^ and Q is quite high, i.e. near the disk midplane. For 
these reasons, we henceforth consider a late stage of solar 
nebula evolution in which the gaseous disk has depleted by 
a factor of 10 (/s — 10~^) and settling of dust grains onto 
the midplane has proceeded signiflcantly {fd < 10~^). 

Let us first estimate the thickness of this dust-rich layer. 
We assume that C has increased from the well-mixed value 
of 10~^ to 50. Taking into accoimt the reduction in gas 
density, the vertical scale height of dust grains is expected 
to decrease by a factor of 500 compared to the gas scale 
height of 0.1 AU at ii = 2 AU. Thus a thickness of 3 x 10* 
km is estimated. This thin dust-rich layer around the mid- 
plane will be mostly magnetically dead because of the ex- 
tremely high dust density there, and the resulting high 
recombination rate on grain surfaces (Sano et al. 2000). 
Yet, the thin surfaces of this dust layer will be sufficiently 
ionized by cosmic rays and stellar X-ray radiation to sup- 
port the MRI and generate current sheets, considering that 
there are now few dust grains suspended in the upper at- 
mosphere. Current sheets that form in these thin layers 
will ionize nearby regions within the dust lay(T, into optical 
depth of order unity. Dust grains are heated after passing 
through these current sheets, sometimes up to their melt- 
ing temperature (Tmeit ~ 1600 K) as shown in § 3.1, and 
then radiate strongly in the near-infrared, subsequently 
raising the temperature of dust grains in the surrounding 
environment. Collisions with hot dust grains may then 
ionize* the gas in an initially magnetically dead area above 
the critical value, Xcrit- 

However, the integrated optical depth of the protosolar 
nebula at i? = 2 AU to the near-infrared radiation in the 
vertical direction is estimated to be r = «;/_rS ~ 2 x 10^ 
(D'Alessio, Calvet, & Hartmann 2001; Chiang et al. 2001). 
Even if we assume that some of the smallest, submicron- 
sized grains have been removed from the nebula with the 
gas, the magnetically active region aflected by current 
sheets would still be only a few percent of the dust-rich 
layer. To melt the majority of the dust grains that reside 
in the midplane and form chondrulcs, the entire layer must 
be in motion so that all dust particles are exposed in the 
magnetically active surface layer at some point. 

When most of the solids have settled to the midplane, 
vertical shear between the gas and the differentially ro- 
tating, dense particle layer can drive turbulence within 
a 10*- 10"'' km thick boundary layer (Cuzzi, Dobrovolskis, 
& Champncy 1993; Cuzzi, Dobrovolskis, & Hogan 1996). 
Although vertical velocity gradients vanish by symmetry 
very near the midplane, a turbulence model that self- 
consistently deals with generation, transport, and damp- 
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inside current sheets is given by (Chiang et al. 2001) 



Stt^X / 2Bx{T^eit)RdR 

JlAV 

'AMexp(9/A^)-l]^"' 



2 X 103 



ergs 



-1 



(19) 



for Tig = 10^^ cm~3 at ii = 3AU. Here, is the wave- 
length of radiation in microns. Taking Ics ~ 10^ km and 
Lbh ~ 0.03 AU, we obtain Ics/Lbh ~ 3 x 10^''' for the 
fraction covered by current sheets in the region between 1 
and 2 AU, the site of their formation in an evolved disk 
thinned by /e = 10""^ (see Figure lb). This crude estimate 
for the excess near-lR emission from current sheets may 
partly explain the observed excesses in the range of 10'^^- 
10^3 erg s~^. Uncertainties in modeling of protoplanetary 
disks at present are admittedly large. Our mechanism in- 
volving current sheets could be taken as an example of 
producing near-lR excesses when magnctohydrodynami- 
cal effects are considered. 

Ultimately, these ideas can be tested by searching 
for distinguishable observable signatures using techniques 
that have recently emerged. Interferometric imaging of 
LkHa in the H and K bands by Tuthill, Monnier, & Danchi 
(2001) has shown a central gap (or cavity) and a hot in- 
ner edge, observed at 3.4 AU from the central star when a 
distance of 160 pc to LkHa is assumed. Earlier theoretical 
predictions for its radius to fit the near-IR inflections in 
the SED (~ lOE*; Hillenbrand et al. 1992) are too small 
by around an order of magnitude. If we can divide the disk 
into several annular regions and obtain photometric mea- 
surements of each annulus, we suggest that there would be 
a local maximum of near-IR radiation between 1 and 2 AU, 
for systems like our primordial solar nebula (although for 
other young stellar objects, the exact position will depend 
on stellar and disk parameters and should be recalculated) . 
In standard protoplanetary disk models such as CG97, the 
mentioned region reaches temperatures of only 100—400 K 
and is expected to emit poorly in the near-IR. In the near 
future, a multi-band, adaptive optics coronagraph under 
development at the American Museum of Natural History 
(Oppenheimer 2000) will be used to directly view proto- 
planetary disks and will present another possibility to test 
these theories. 



current sheet would survive in that environment. To study 
these effects in detail a dynamical model will be necessary. 
This requires modeling the interplay between at least four 
physical processes: pressure imbalance, ionization evolu- 
tion, temperature evolution, and steepening of the mag- 
netic field. These four processes are inherently connected 
and feed back on one another. That means they must 
be solved simultaneously in a self-consistent fashion. Ra- 
diative transfer and a power-law size distribution of dust 
grains must be incorporated in a subsequent step to un- 
derstand chondrule formation in these dynamical models. 
The large range of timescales and length scales involved 
presents another challenge for modeling of this compli- 
cated set of processes. 
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for pointing out the importance of high oxygen fugacities 
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knowledges support by NSF Career grant AST99-85392, 
NSF grant AST03-07793, and NASA ATP grant NAG5- 
10103. D. S. E. acknowledges support from NASA Cos- 
mochemistry grant NAGS- 12855. 



4. CAVEATS 

Our conclusions suS^er from several caveats. All models 
of the current sheets that have been done to date have as- 
sumed a steady state flow (BZ94; Brandenburg & Zweibel 
1995; Zweibel & Brandenburg 1997; Heitsch & Zweibel 
2003). However, the heating and ionization that occur in 
the current sheet are not accounted for self-consistently 
in these models. Both will increase the ion pressure in 
the center of the sheet, cither saturating the singularity 
or even blowing the sheet apart. Also, we invoke turbu- 
lence induced by the MRI when calculating magnetic field 
fluctuations but ignore its effect on the evolution of cur- 
rent sheets, e.g. the converging motion of the neutrals in 
ambipolar diffusion. In our simplified chondrule forma- 
tion model, we let gas and dust travel through a current 
sheet at ~ 1 km s~^, but it is not clear how well this ap- 
proximates the actual condition near a current sheet or if a 
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Fig. 1. — Relevant length scales, (a) Assumed parameter values are = 1, B = 1 G, and xi = lO'^^. The critical wavelength of the MRI, 
resistive diffusion scale, and (twice) the scale height of the disk arc shown as solid, dotted, and dashed lines, respectively. See text for their 
definitions. The shaded region satisfying Af; < \bh < 2/i, where the MRI occurs, extends from 2 AU to 4 AU. (b) Same as (a) but using 
/e = 10~^. This is appropriate for a late stage of the solar nebula evolution. Note that the magnetically active region moves inward to 1-2 
AU. 



z® > 

X 



null 



\/ \/ \/ v.vW'MAA 

Vi . 



JxB 



Vi 
JxB 



Fig. 2. — Configuration for magnetic fields in the one-dimensional geometry considered in this paper. Current density (J) points in the +z 
direction, out of the page. The directions of various terms in Eqns. (3) and (4) are also shown. A current sheet forms along the magnetic 
null, shown as a dotted vertical line along the center of the figure. 



Chondrule Formation by Current Sheets 



11 



0.004 



1 — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — r 




-2000 



2000 



O 
O 



-6 - 



101 



IQi 



_L 



2 
LOG X (km) 



Fig. 3. — (a) Heating profile for our fiducial model with /s = 0.1. Note the symmetry about a; = 0. Ohmic dissipation (dotted line) is 
concentrated near center of current sheet, while the ambipolar diffusion heating (dashed line) is centered off-peak. The latter exceeds the 
contribution from the former by a factor of ~6 in our fiducial model, (b) Heating profiles as a function of Bo in a log-log scale. In this 
case, the gas density (rig) was varied, and each curve is labeled with the corresponding value of Ug used. The peak heating rate increases 
approximately inversely proportional to the square of gas density; the width of current sheet varies as described by Eqn. (13). 
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Fig. 4. — Time evolution of gas, dust, and chondrule temperatures at a fixed position in space. Common parameters used for the two runs 
are: Ug = 10^^ cm~^, = 1 /im, Oc = 1 mm, and Bmax = 3 G. Solid lines are used for the run with C = 5 and dotted lines for the run 
with C, = 50. Only for the higher value of C, do dust grains and chondrules reach the melting temperature. We stop running the code when 
the particles evaporate. The initially steep rise in Tg^s is due to intial conditions that are not set self-consistently. 



